Role of three-body interactions in formation of bulk viscosity in liquid 
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With the aim of locating the origin of discrepancy between experimental and computer simulation results on 
bulk viscosity of liquid argon, a molecular dynamic simulation of argon interacting via ab initio pair potential 
and triple-dipole three-body potential has been undertaken. Bulk viscosity, obtained using Green-Kubo 
formula, is different from the values obtained from modeling argon using Lennard- Jones potential, the former 
being closer to the experimental data. The conclusion is made that many-body inter-atomic interaction plays 
a significant role in formation of bulk viscosity. 
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I. INTRODUCTION 

Argon above its melting temperature is a typical sim- 
ple fluid. Consisting of spherical atoms that interact via 
short-range repulsion and long-range attraction, and are 
heavy enough for the quantum effects to be small, fluid 
argon and heavier noble gases are an excellent choice of 
a real system to be used for testing various approaches 
in classical theory of fluids. 

An inter-particle interaction in argon is commonly 
represented by a well known 12-6 Lennard- Jones pair 
potential 1 , 



«lj (r) = 4e L j 



(?)-(? 



in 



The two parameters, ctlj and are usually determined 
by fitting thermodynamic properties, derived from the 
potential (JTJ) by theoretical or computational methods, 
to corresponding experimental data. 

It is known that Lennard- Jones potential is only an 
approximation to real interaction in argon. Several ex- 
perimental results obtained for argon at large pressures 
are better explained if a larger steepness, compared 
to Lennard- Jones, of argon-argon interaction potential 
at small inter-atomic separation distances is taken into 
account 2,3 . Accurate argon-argon interatomic potentials 
have been calculated by direct ab initio quantum chemi- 
cal calculations^ - — or obtained by inversion of experimen- 
tal data 7 -. Moreover, many-body dispersion, exchange 
and induced polarization contributions to inter-atomic 
interactions are not small and noticeably influence ther- 
modynamic properties of argon^ 9 -. The most widely used 
of these contributions is triple-dipole dispersion interac- 
tion, derived by Axilrod and Telle r 10 i 11 and Mutoi 2 -, and 
account of this contribution in addition to ab initio pair 
potential is sufficient to describe thermodynamic proper- 
ties of argon with good accuracy-i 3 - - — . 

By virtue of Henderson theore m 18 ' 19 , which states 
that, for fluids with only pairwise interactions, and under 
given conditions of temperature and density, the pair po- 
tential which gives rise to a given radial distribution func- 



tion g(r) is unique up to a constant, the thermodynamic 
properties of the system with many-body interactions can 
be described by a model system with an appropriate ef- 
fective pair potential. Generally, the effective potential 
depends on the thermodynamic state of the system and 
thermodynamic property to be described 20-22 . Van der 
Hoef and Madden 21 have demonstrated that the account 
of triple-dipole and dipole-dipole-quadrupole dispersion 
interactions moves the effective potential of argon to- 
wards Lennard- Jones form (fTJ. Moreover, the possibility 
of consistent description of many thermodynamic prop- 
erties of argon, using Lennard- Jones potential in a wide 
domain of thermodynamic states^ 3 . - — , suggests that the 
state dependence of the effective potential is weak. 

There is no analogous reason for kinetic properties of a 
system with many-body interactions to be equivalent to 
those of a system with a corresponding effective pair po- 
tential. Nevertheless, experimental data on self-diffusion, 
shear viscosity and thermal conductivity coefficients of 
argon have been shown to be accurately described by 
Lennard- Jones model with the parameters obtained by 
fitting thermodynamic dat a 26 ' 27 . 

Bulk viscosity is a noticeable exception. Bulk viscos- 
ity of argon has been measured experimentally^ - —, and 
its behavior can be qualitatively described by the results 
of a molecular dynamics simulation of a Lennard- Jones 
system 3 ^. However, when results of simulations with 
Lennard- Jones potential are rescaled in an attempt to 
describe experimental data liquid argon, bulk viscosity, 
contrary to other kinetic properties, appears strongly un- 
derestimated (e.g. up to 50% in Ref. [271)]. 

In view of the above, I propose that the source of this 
discrepancy may lie in neglect of many-body interactions. 
Previous molecular dynamics simulations of systems con- 
sisting of 108 particles interacting via ab initio pair po- 
tential and Axilrod- Teller-Muto (ATM) interaction indi- 
cated that a triple-dipole interaction does not affect the 
bulk viscosity of liquid xenon near its triple point— and 
dense gaseous krypton^. However, the error in the val- 
ues of bulk viscosity obtained from molecular dynamics 
simulation of the systems with such a small number of 
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particles can be quite large. For example, the values of 
the reduced bulk viscosity of the Lennard- Jones systems 
consisting of 128 and 256 particles at the reduced temper- 
ature T* = 0.722 and the reduced density p* = 0.8442, 
reported in Refsl36l. l39H4ll . range from 0.89 to 1.47, with 
the ratio of the latter to the former of 1.65. 

This paper presents the results of more accurate molec- 
ular dynamics simulations of a liquid consisting of 1372 
argon atoms with ab initio+ATM interaction, which 
demonstrate that bulk viscosity, determined from Green- 
Kubo formulae, significantly changes with the account 
of three-body interaction, moving results towards exper- 
imental data. 



II. INTERACTION 

Nasrabad et al 16 undertook a Monte Carlo simulation 
of argon using combination of ab initio pair interaction 4 
and ATM triple-dipole dispersion interaction 1 ^ to test 
their ability to predict vapor-liquid equilibrium. Al- 
though more accurate ab initio pair potentials for ar- 
gon have become available recently^, and other many- 
body contributions to inter-atom interaction can be 
calculated^, we use the same interaction as Nasrabad et 
al because, being able to predict accurately the phase 
diagram of argon 16 , it is computationally more efficient. 

Specifically, the ab initio pair interaction potential 
used in the present work is described by a function^ 6 - 

u 2 (r)=Ae-«r +^ 2 + ^/ 2nM )^, (2) 

n — 3 

where 

/ 2 „M) = i- e - b ''£^, (3) 

and numerical values of the parameters A, a, ft, b, and 
C 2n are given in Ref. [16j. The ATM triple-dipole inter- 
action has formiS 

, . 1 + 3 cos a cos /3 cos 7 . . 
u 3 (r 12 ,r 23 ,r 31 ) = v , (4) 

'12'23'31 

where the are the lengths of the sides, a, j3, and 7 are 
the angles of the triangle formed by three argon atoms, 
and v = 7.32 • 10~ 108 J-m 9 for argoni 3 -! 4 .. 

For simulations of argon using Lennard- Jones potential 
IHJ the values ct lj = 3.3952 A and e LJ = 116.79 K are 
used 2 ^. 



III. SIMULATION 

Meier et a' 36 undertook a systematic study of the in- 
fluence of the number of particles and the cutoff radius 
for pair interaction on the bulk viscosity of Lennard- 
Jones system. In view of their results, simulations were 



performed in a cubic box containing N = 1372 parti- 
cles, and the cutoff radius for pair interactions was set 
to 5<7lj. Three-body interactions were cut off when the 
distance between any pair of the atoms in the triplet ex- 
ceeded one quarter of the simulation box length (around 
3ol,t for the densities studied in this work). Usual pe- 
riodic boundary conditions and minimum image conven- 
tion were applied. The simulations were started with the 
particles in a face-centered-cubic lattice, with randomly 
assigned velocities. Forces arising from three-body in- 
teractions were calculated using formulas given by Allen 
and Tildesley 4 ^, and an expression for forces due to ab 
initio pair interaction was obtained by applying gradient 
operator to Eq. ([2]). Newton's equations of motion were 
solved using velocity- Verlet algorithm with the time step 
At • v^j7to/«7Lj = 0.003. 

The runs were made at the experimental densities at 
various temperatures along the 40 atm isochore, taken 
from Rcf. [33]. Every simulation was initiated in the 
NVT ensemble and run for at least 2-10 5 time steps to 
attain thermodynamic equilibrium. After equilibration 
the thermostat was turned off and the NVE ensemble 
was invoked to calculate bulk and shear viscosities. The 
length of the production period was 4-10 6 time steps for 
the system interacting via Lennard- Jones potential, and 
between 10 6 and 3T0 6 time steps for the system with ab 
initio + ATM interaction, depending on the state point. 

Bulk viscosity, £, and shear viscosity, 77, were calculated 
using Green-Kubo formulas^ 3 -: 

f = rW° (MWo)) dt, (5) 



where V is volume, ks is Boltzmann constant, T is tem- 
perature, t is time, Sp = p — (p) is the deviation of the 
instantaneous pressure p from its average value (p) , cr Q( 3 
is an off-diagonal element of the stress tensor, the angu- 
lar brackets denote equilibrium ensemble averages over 
short trajectory sections of the phase-space trajectory of 
the system with multiple (every time step) time origins 
to . The stress tensor was calculated using formulae given 
by Lee and Cummings The integration in Eqs ([S]) and 
([6]) was carried out up to tl = LJ c, where L is simulation 
box length and c is sound velocity taken from Ref. (33| . 
Depending on the state point, the value of t£ was be- 
tween 4.80 and 11.25 ps. The statistical error in time 
correlation functions was estimated using formula given 
by Frenkel and Smit 4 ^, 

a({X(t)X(0)))^ x f^(X"(0)), (7) 

V ''run 

where i run is the length of the simulation, and the cor- 
relation time Tx was approximated as the time during 
which time correlation function decays e ~ 2.718 times. 
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Figure 1. Bulk viscosity of liquid argon at T = (90—140) K. 
Error bars connected with solid and dashed lines correspond 
to the simulation results with ab initio + ATM and Lennard- 
Jones interaction, respectively. Experimental points are taken 
from Refs [33| (circles, pressure 40 atm) and [29| (square with 
error bar, pressure 40 kg/cm 2 ). 



Figure 3. Shear viscosity of liquid argon at T — (90—140) K. 
Error bars connected with solid and dashed lines correspond 
to the simulation results with ab initio + ATM and Lennard- 
Jones interaction, respectively. Dotted line corresponds to the 
interpolation data for pressure 40 atm taken from Ref. [Hol |. 
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Figure 2. Time-correlation functions C(t) used for calcu- 
lation of bulk viscosity at density 1.258 g/cm 3 . Solid and 
dashed lines correspond to the simulation results with ab ini- 
tio + ATM and Lennard- Jones interaction, respectively. 



IV. RESULTS 

Fig. [T] and Table U present simulation results for the 
bulk viscosity obtained using ab initio + ATM (Eqs @ 
and (QJ) and Lennard- Jones (Eq. ([1])) interaction, re- 
spectively. Bulk viscosity, determined from Green-Kubo 
formulas, changes with the account of three-body in- 
teraction, moving towards experimental data. How- 
ever, this change is not sufficient to obtain numerical 
agreement with experiment, especially at lower densities. 
Typical behavior of time correlation functions C(t) = 
(Sp(t)5p(0)) is shown in Fig. 03 



Fernandez et a? 7 demonstrated that, contrary to bulk 
viscosity, the values of shear viscosity of argon ob- 
tained from molecular dynamics simulation of a Lennard- 
Jones system agree with experimental data. Lee and 
Cummings 44 and Marcelli et al 4 ® found that the influence 
of triple-dipole interaction on shear viscosity of argon is 
small. The results of the present simulation, shown in 
Fig. 02 and Table [Q agree with these findings. 



V. CONCLUSION 

The message of this paper is that many-body inter- 
actions play a more substantial role in determining the 
value of the bulk viscosity than other transport coeffi- 
cients. The present results from the molecular dynamic 
simulation of liquid argon demonstrate that even account 
of a single many-body contribution, ATM triple-dipole 
interaction, shifts the values of the bulk viscosity of argon 
towards experimental data. Larger sensitivity of the bulk 
viscosity to many-body interaction, compared to other 
transport coefficients, can be intuitively explained in the 
case of gaseous state. Bulk viscosity of a non-relativistic 
monoatomic gas calculated from the Boltzmann equa- 
tion, which takes into account only pair collisions of 
atoms, appears to be zero, in contrast to heat conductiv- 
ity and shear viscosity which have non-zero values in the 
same approximation 47 . A non-zero value of bulk viscosity 
appears in the approximations corresponding to higher- 
order terms in the virial expansio n 48 ! 49 , which correspond 
to the explicit account of at least three-atom collisions 
which, in turn, are sensitive to three-body inter-atomic 
interaction. 
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T K 


n P'/ptti^ 
Hi &/ , - ,±li 


Bulk viscosity f, mps 


Shear viscosity r\, mps 


LJ 


AI±ATM 


Ref. [33] 


LJ 


AI±ATM 


Ref. [33] 


Ref. [50] 


90 


1.390 


1.10 ±0.04 


1.54 ±0.10 


1.82 


2.31 ±0.04 


2.44 ± 0.07 


2.33 


2.57 


100 


1.327 


1.03 ±0.03 


1.48 ±0.09 


1.57 


1.78 ±0.03 


1.87 ±0.06 


1.86 


1.92 


110 


1.258 


1.04 ±0.02 


1.33 ±0.05 


1.39 


1.38 ±0.02 


1.39 ±0.03 


1.51 


1.48 


120 


1.182 


0.99 ±0.03 


1.35 ±0.05 


1.51 


1.09 ±0.02 


1.12 ±0.02 


1.19 


1.15 


130 


1.092 


1.04 ± 0.04 


1.32 ±0.10 


1.71 


0.86 ±0.02 


0.87 ±0.03 


0.88 


0.89 


135 


1.037 


1.03 ±0.04 


1.21 ±0.10 


1.93 


0.73 ± 0.02 


0.70 ± 0.03 


0.760 


0.77 


140 


0.968 


1.00 ±0.05 


1.12 ±0.10 


2.53 


0.65 ±0.02 


0.65 ±0.03 


0.642 


0.65 



Table I. Bulk and shear viscosities of argon obrained from molecular dynamics simulations using Lennard- Jones (LJ) and ab 
initio pair + Axilrod-Teller-Muto three-body (AI±ATM) interaction, and corresponding experimental dat ai 33 ' 50 . Error in the 
simulation data is calculated using Eq. ([7}. 
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